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Abstract 

Motivated by recent observations that suggest a low density of old stars around the Milky Way supermassive 
black hole, models for the nuclear star cluster are considered that have not yet reached a steady state under 
the influence of gravitational encounters. A core of initial radius 1 — 1.5 pc evolves to a size of approximately 
0.5 pc after 10 Gyr, roughly the size of the observed core. The absence of a Bahcall-Wolf cusp is naturally 
explained in these models, without the need for fine-tuning or implausible initial conditions. In the absence of 
a cusp, the time for a IOMq black hole to spiral in to the Galactic center from an initial distance of 5 pc can be 
much greater than 10 Gyr Assuming that the stellar black holes had the same phase-space distribution initially 
as the stars, their density after 5 — 10 Gyr is predicted to rise very steeply going into the stellar core, but could 
remain substantially below the densities inferred from steady-state models that include a steep density cusp 
in the stars. Possible mechanisms for the creation of the parsec-scale initial core include destruction of stars 
on centrophilic orbits in a pre-existing triaxial nucleus, inhibited star formation near the supermassive black 
hole, or ejection of stars by a massive binary. The implications of these models are discussed for the rates of 
gravitational-wave inspiral events, as well as other physical processes that depend on a high density of stars or 
stellar mass black holes near SgrA*. 



1. INTRODUCTION 

Near-infrared imaging reveals a cluster of stars around the 
~ 4 X 1O^M0 supermassive black hole (SMBH) at the cen- 
ter of the Galaxy. Much recent work has focussed on the 
youn g stars that dominate the total hght in the central par - 
sec (F orrest et all 19871: 1 Allen et alJ199alKrabbe et al.ll991h . 
These stars have masses of 10 — 6OM0 and appear to have 
forme d in one or m ore starbursts during the last few million 
years (i Krabbe etal] ri995: Paumard et al. 2006). While their 
total numbers are small, the density of the young stars in- 
creases very steeply toward the SMBH jGenzel et al J 120031 
iPaumard et alj|2006l) . How such massive stars can form so 
deep in the SMBH potential well remains an unsolved prob- 
lem (pilexander 2005; Paumai'd 2008). 

The dominant population at the Galactic center consists 
of old stars: mostly metal-rich, M, K a nd G type giants 
with masses of one to a few Solar masses jBlum et 
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iDavies et al.ll2009l) . The late type giants dominate the total 
flux outside the central parsec, and seeing-limited observa- 
tions as early as the 1960s showed that their surface den- 
sity follows a power law, E ^ R^^ at R <^ 10 pc, imply- 
ing a space density profile » ^ jBecklin & Neugebauej 
mm iMcGinn et al.lll989l: iHaller et al.lll996^ Inside ^ 0.5 
pc, a drop in the CO absorption strength dSellgren et alj 
Il990t iHaller et alj |1996) sign als a decrea s e in the projected 
density of old star s (Genzel et al.lll996t iFiger et alj [20001: 
IScoville et all 120031: [Figer et al.l l2003h . However the possi- 
bility of contamination by light from early type stars at these 
radii made such inferences uncertain. 

Three recen t studies have clarified the situation. 
iBuchholz et al.l (|2009l) took deep narrow band images 
of the inner parsec and used a CO absorption feature to 
distinguish early from late type stars. They classified roughly 
3000 stars down to a /T-band magnitude of ~ 15.5. Of 
these roughly 300 were early types, distributed as a steep 
power-law around the SMBH. The late type stars showed a 
very different distribution: a core of radius ^ 0.5 pc, with a 



constant or even declining surface density inside. Do et al. 
(2009) obtained higher dispersion spectra of a smaller sample 
of stars in a set of fields within 0.16 pc from SgrA*, estimated 
to be 40% complete down to niK = 15.5. They also found 
a flat or centrally-declining projected density of old stars. 
iBartko et al. (2009b) found a flat distribution of late-type 
stars with nik < 15.5 based on spectroscopic identifications in 
a set of fields in the central 25". 

Inferring the space density profile of the old stars from the 
number counts in the inner parsec is difficult, since the pro- 
jected density is apparently dominated by stars that are far 
from the center But the behavior of the number counts is dif- 
ficult to reconcile with models in which the space density in- 
creases inwards. Buchholz et al. (2009) found a best-fit power 
law of E o<: R'^-^'' (projected density) for late type stars inside 
6" and E Z?^"^ outside 6". This corresponds to a space den- 
sity that drops toward the center, although the counting statis- 
tics were also consistent with a flatter profile. Do et al. (2009) 
inferred a steeper rate of central decline, E 0= 7?"^^, from their 
smaller data set, though with larger uncertainties. Do et al. 
compared the numbers counts with projections of power law 
models in n(r) and concluded that the density could not in- 
crease faster than « ^ ' toward the center, although the pre- 
ferred dependence was shallower. 

As the authors of these three studies note, previous descrip- 
tions of the stellar density as a broken power law, with n ~ 
r~'-^ — r^^-^ inside ^ 0.4p c and n ^ r^^ outside dGenzel et alj 
I200I ISchHdel et alj|2007l) . were somewhat misleading, since 
the counts at small radii were dominated by early type stars 
while the counts at large radii were dominated by late type 
stars. Presumably, the physical mechanisms responsible for 
creating these two populations, and placing them on their cur- 
rent orbits, were different and occurred at different times. 

Models with a central "hole" naturally resolve one long- 
standing puzzle: why virial estimates of the SMBH mass 
based on velocitie s of late type stars ( Ghez et al . 1 9 91; 
iGenzel et all l2000l iChakrabartv & Sahal 126011: lEckartetal] 
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I2OO2I) gave systematically lo wer values than the mass inferred 
from the inner S-star orbits dGillessen et alj|2009t iGhez et alj 
I2OO8I) . If most of the old stars within the projected central 
parsec are actually far from the SMBH, their motions will be 
relatively unaffected by its gravitational force, yielding spuri- 
ously low virial masses. Recent studies that explic itly incor- 
porate a low central densi ty for the late type stars jZhu et alj 
I2008t ISchodel et al.ll2009l) yield virial masses for the SMBH 
that for the first time are consistent with the mass obtained 
from the orbital fits. 

However the low density of old stars at the Galactic cen- 
ter creates new puzzles. A standard assumption has long 
been that stars near the SMBH should exhibit a cusp in the 
density, n - r^^, y w 7/4. This is the Bahcall-Wolf (1976) 
solution, which holds for a relaxed (in the sense of gravi- 
tational encounters) population of stars moving in a point- 
mass potential. A number of theoretical studies have argued 
that the Bahcall-Wolf solution is a robust outcome, depend- 
ing only weakly on the initial conditions and the range of 
mass groups present jMurphv et al.|[T99Tt iFreitag et al.ll2006l : 
iHopman & Alexandeij|2006h . 

It is possible that a Bahcall-Wolf cusp is present but that 
the stellar luminosity function changes insid^ 0-5 PC. such 
that the brightest st ars are missin g (iBuchholz et alJ 120091 : 
iDo et all 12009: Bartko et al.ll2009bl) . This could be a result 
of (physical) collisions before or during the red giant phase, 
which strip stellar env elopes and keep stars from reaching 
their peak luminosities (iGenzel et al. 1996: Alexander 1999: 
iBailey & D avies"1999^: or tidal interactions between single 
stars and the SMBH ( Davies & Kin g 2005): or it could result 
from an initial mass function that is truncated below ^ SMp,, 
since these are the stars that w ould otherwise domi nate the K- 
band number counts now (e.g.l Navakshin & Sunv aev 2005). 

While a "hidden" density cusp is possible, it seems reason- 
able also to consider models in which the observed stars are 
representative of the unobserved stars. In these models, both 
populations would have a relatively low density in the inner 
parsec. Aside from Occam's principle, a number of other mo- 
tivations exist for considering such models: 

1 . The relaxation time at the Galactic center is long. The 
stellar mass density at a distance of 2 — 3 pc from SgrA*- 
roughly the SMBH influence radius rinfl - implies a two-body 
relaxation time (for Solar-mass stars) of 20-30 Gyr This is 
the relevant time scale for refilling of an evacuated core if its 
initial radius is ~ nnfl. 

2. Stellar kinematics suggest a low mass density in the inner 
parsec. Schodel et al. (2009) found that proper motion data 
were consistent with a range of models for the distribution of 
mass in the inner parsec, but the best-fitting models had a flat 
or declining mass density toward SgrA*, similar to what is 
seen in the number counts of the old giants. 

3. Physical collisions fail, by a wide margin, to predict 
the observed depletion of giant stars in the faintest magnitude 
bins. Dale et al. (2009) concluded that giants in the 15 > 
niK > 12 magnitude range would only be signficantly depleted 
within ^ 0.01 pc, even assuming a density of unseen colliders 
(main sequence stars, stellar-mass BHs) that was four times 
larger than in the dynamically relaxed models with a Bahcall- 
Wolf cusp (e.g. Hopman & Alexander 2006). 

4. Whatever their origin, cores are ubiquitous components 
of galaxies with SMBHs, at least in galaxies that are bright 
enough or near enoug h for parsec-scale features to be resolved 
dFerrarese et al.l2 006). Core radii are roughly equal to SMBH 
influence radii (e.g. ,Graham.2004.) . which in the case of the 



Milky Way would predict a core of radius 2 — 3 pc. 

This paper examines the viability of models for the Galac- 
tic center that include a low-density core. Core models are 
examined first from a structural point of view (§2), then from 
the point of view of self-consistent equilibria (§3), and finally 
from an evolutionary standpoint (§4). The basic result is that 
core models "work": they reproduce the number counts and 
kinematics of the late type stars, without invoking physical 
collisions (collisions would be extremely rare in these mod- 
els) and without the necessity of fine-tuning. We show that a 
core of the size currently observed is a natural consequence of 
two-body relaxation acting over 10 Gyr, starting from a core 
of radius ^1 — 1.5 pc; the n ^ r^^ density profile outside the 
core gradually extends inward as the core shrinks and as the 
stellar density evolves toward, but does not fully reach, the 
Bahcall-Wolf form after 10 Gyr. 

The relaxation time that sets the rate of evolution in these 
models depends inversely on the mean stellar mass, roughly 
IMq under standard a ssumptions about the initial mass func- 
tion (lAlexandeiil2005) . But an old population also contains 
stellar remnants, including ^ lOM© black holes. In the ab- 
sence of a core, black holes initially at distances of 4 — 5 pc 
from SgrA* would spiral all the way in to the center after 10 
Gyr (lMorrisll993h . potentially dominating the total mass den- 
sity inside ~ 10^^ parsec. In the models considered here, the 
black hole orbits tend to decay no farther than the core ra- 
dius as determined by the stars (§5). The result, after sev- 
eral Gyr, is a rather different distribution of black holes than 
in the coUisionally-relaxed models. This difference has po- 
tentially important implications for the rates of gravitational 
wave driven inspirals or for other dynamical processes that 
postulate a dense cluster of black holes around SgrA*, as dis- 
cussed in §6. 

The observations motivating this paper are fairly new, and 
if recent history is any guide, our observational understanding 
of the Galactic center will continue to change. It is conceiv- 
able that the low apparent density of late type stars is an ar- 
tifact due to faulty stellar classifications, improper treatment 
of extinction, confusion, or some other factor Stars slightly 
fainter than the current limit for robust detection (mK ~ 15.5) 
may turn out, once detected, to have a very different spatial 
distribution, more similar to that predicted by Bahcall and 
Wolf (1976). Our theoretical understanding of stellar colli- 
sions may also change; for instance, new, more effective chan- 
nels for removal of stellar envelopes could be discovered. Any 
of these developments would lessen the relevance of the mod- 
els discussed here to the center of the Milky Way. However, 
core models are also applicable to other galaxies that con- 
tain SMBHs, many of which have cores that can not plausibly 
be explained except in terms of a general depletion of stars. 
And a robust conclusion to be drawn from this work is that 
the distribution of stars and stellar remnants at the center of 
the Milky Way should still reflect to some extent the details 
of the Galaxy's formation - though the imprint of the initial 
conditions may turn out to be less extreme than in the models 
considered here. 

We assume throughout a distance to SgrA * of 8.0 kpc 
(lEisenhauer et al.ll2005t iGillessen et ani2009l) . At this dis- 
tance, one arc second corresponds to a linear distance of 
0.0388 pc and one parsec corresponds to 25.78". We 
also fix the SMBH mass to 4.0 x IO^Mq (iGhez et al J 120081: 
IGillesse n et al. 2009). 

2. PROPERTffiS OF THE STELLAR DISTRIBUTION 
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2.1. Number density 

As summarized above, a relatively low density of late-type 
stars near the center of the Milky Way was independently in- 
ferred by Buchholz et al. (2009), Do et al. (2009) and Bartko 
et al. (2009). The first of these three studies was based on 
the largest sample of stars, and we adopt the Buchholz et al. 
number counts as basis for the discussion that follows. 

Buchholz et al. (2009; hereafter BSE09) took deep, narrow- 
band images {K < 15.5) in the near-IR K band and distin- 
guished late type (old) from early type (young) stars using CO 
equivalent widths, calibrated from existing samples. Extinc- 
tion was estimated star by star by comparing observed SEDs 
with a blackbody of variable extinction. BSE09 identified a 
"quality 1" sample consisting of late type stars in which the 
CO band depth exceeded the cutoff value for early type stars 
by lo or more, and which did not fall into any of the other 
classes defined by them (AGB stars, foreground sources, very 
red objects). The resulting data set contains 2955 stars down 
to a magnitude limit of K — 15.5. Figure [T] shows the BSE09 
number counts (their Fig. 11) plotted versus projected dis- 
tance from SgrA*. 

The BSE09 data are believed complete to a limiting mag- 
nitude K = 15.5 within a projected distance of ^ 20" from 
SgrA*. Beyond this radius, the early type stars are insignifi- 
cant in numbers compared with the late type stars. Schodel et 
al. (2007) present total number counts in a more extended 
region, derived from the ISAAC NIR camera on the VLT. 
These counts extend to ^ 60" and to a limiting magnitude 
of niK w 17. Down to the magnitude limit (K = 15.5) of the 
BSE09 sample, the ISAAC counts are ^ 80% complete be- 
yond - 15". 

Figure [T] includes number counts from all stars in the 
ISAAC sample with K < 15.5, and within the region 20" < 
R < 60". The observed positions were converted to projected 
densities using the kernel routine described in Schodel et al. 
(2007) including corrections for completeness and crowding, 
as discussed in that paper. Confidence intervals were derived 
via the bootstrap. The ISAAC counts are seen to match very 
well onto the BSE09 counts, without the necessity for any ad 
hoc adjustment of the normalizations. 

At larger radii, out to R ^ lOpc 200", a number of 
studies have found n ^ r^^-^ for the distribution of stel- 
lar light fee cklin & Neugebauer 1968; Halleretal. 1996; 
iGenzelet all [20001: ISchodel et al.ll2007l) . In all of what fol- 
lows, we assume an asymptotic slope dlogn/dlogr = —1.8. 

This leads us to parametrize the space density of the old 
stellar population as a power law with an inner core: 
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with Y= 1.8. The parameter a controls the sharpness of the 
transition from outer to inner slopes. This model was fit to the 
number counts after projecting it onto the plane of the sky; in 
this way one avoids fitting a model to the surface density that 
does not correspond to a physical space density. 

As is often found when modelli ng the lumin osity profiles of 
cored elliptical galaxies (e.g. Ter zic & Grah am 2005), large 
values of a gave the best fits. The thick curve in Figure [T| 
shows the best-fit model when a = 4; it has 

no==0.21pc--\ ro = 5.5" = 0.21 pc, Y, = -1.0. (2) 

The negative value of y, implies a space density that decreases 
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Fig. 1. — Upper panel: azimuthally-averaged surface density of late type 
stars at the Galactic center. Open circles are number counts of the "qualtity 1" 
late type stars brighter than = 15.5 in the sample of Buchholz et al. (2009) 
(their Fig. 11). Filled circles are derived from a kernel estimate of the pro- 
jected density of all stars with K< 15.5 and R > 20" in the sample of Schodel 
et al. (2007), after connections for crowding and completeness. Curves show 
fits of the projected, parametric model, equation (T), to the number count 
data. Heavy curve is the best fit when the inner power-law slope is uncon- 
strained; dashed curve shows the best fit when the inner slope is set at —0.5, 
the shallowest profile consistent with an isotropic velocity distribution. Both 
fits assume a = 4 and y = 1.8. Lower panel: space density profiles corre- 
sponding to the projected profiles in the upper panel. 



toward the center (lower panel of Fig. [TJ - a central "hole." 

Fits of this model to the data were found never to be ter- 
ribly good > 17). There is a broad local maximum in 
the number counts at /? « 15" which would require additional 
parameters to fit. (In the 2d counts. Figure 13 of BSE09, this 
overdensity is seen to be roughly symmetric about the origin.) 
The central minimum in the projected density is also difficult 
to reproduce. 

Figure |2] shows the distribution of inner slopes y, obtained 
from the fits to lO'* Monte-Carlo samples bootstrapped from 
the number count data in Figure [T] 90% of the values fall in 
the range 

- 3.5 <y, < 0.82. (3) 
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defines the local relaxation time as 
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Fig. 2. — Distribution of inner slopes from fits to 10"* Monte-Carlo samples 
bootstrapped from the number count data in Fig.[T] Dashed lines delineate a 
90% confidence interval. 



Negative values of y, (corresponding to centrally-decreasing 
densities) are clearly preferred, though positive values are also 
acceptable. 

If the stellar velocity distribution is isotropic, the space den- 
sity must increase at least as fast as r^"-^ near the SMBH (§3). 
The dotted curve in Figure[T|shows the best fit when y, is fixed 
at 1/2. The other parameters are 



no = 0.12pc-^ ro = 7.7" = 0.30 pc. 



(4) 



This model can be (crudely) thought of having the minimum 
central density consistent with an isotropic velocity distribu- 
tion. This fit is only slightly worse in a sense (19.0 vs. 
17.5) than the fit with unconstrained inner slope, although it 
implies a very different space density profile inside ~ 0.5 pc 
(Fig- [II lower panel). 

A standard definition of the core radius is the projected ra- 
dius where the surface density falls to 1 /2 of its central value 
(e.g. King 1962). As a practical definition of the "central 
density," we take the value at 1" w 0.04 pc projected radius. 
Based on this definition, the core radius of the unconstrained 
fit in Figure [1] is 



15.2" w 0.59 pc 



(5) 



and the core radius of the constrained fit (y, = —1/2) is 

rcore ~ 10.8" w 0.42 pc. (6) 

Thus, the core radius of the Milky Way nuclear star cluster is 
- 0.5 pc. 

2.2. Relaxation time 

The relatively low central density of stars implied by Fig- 
ure [T] and the apparently non-relaxed form of n{r), are sug- 
gestive of a long two-body relaxation time. Spitzer (1987) 
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where o is the rms velocity in any direction, p is the stel- 
lar mass density, m is the mass of one star, and In A is the 
Coulomb logarithm. This expression assumes that all stars 
have the same mass and that their velocity distribution is 
isotropic (Maxwellian). 

To apply equation (O we need an estimate of p(r). Dynam- 
ical estimates of p are to be preferred, given the large system- 
atic uncertainties associated with converting a l uminosity den- 
sity into amass density near the Galactic center (ISchodel et alJ 
l2007tlBuchholz et alJl2009 ). 

Inside ^ 1 pc, the best dynamical constraints on p come 
from the recent proper motion study of Schodel et al. (2009). 
These authors detected, for the first time, an unambiguous sig- 
nature of the gravitational force from the distributed mass on 
the stellar motions, in the region 0.25 ^ r ^ 1 pc; inside this 
region the gravitational force from the SMBH can not be dis- 
entagled from that of the stars. Schodel et al. inferred a dis- 
tributed mass of - 1 ±0.5 X IO^Mq in a sphere of radius 1 
pc around SgrA*. The functional form of p(r) was not well 
constrained, although the formal best fits were obtained with 
models in which the mass density decreased toward the center. 

We use the Schodel et al. results to normalize the density 
of our "isotropic" n(r), equations ([T]|4l). Defining 

M,^(r<lpc) 

M^= ^ r (9) 

the mass density becomes 

p(r) = 9.9 X lO^Mopc-^M,^-''-^ (l +^4)""-^^^ (10) 

with ^ = r/0.30 pc. We then calculate the (isotropic) velocity 
dispersion from the Jeans equation, 

p(r)o(r)2 = G^ a'rV'-2[M.+M4<r')]p(r'). (11) 

Figure [3] shows the resulting f,-(r) assuming InA = 15, m = 
Mq, and three values of M*. 

The Schodel et al. (2009) proper motion data do not con- 
strain the density beyond — 1 pc. In this region, we assume 
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a simple continuation of the model fit at smaller radii. Below 
we argue for a "preferred" density model with p = 1 .5, but in 
Figure[3]we allow the normalization: 



p(lpc) 
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to have the values p = (0.75, 1.5,3). The Jeans equation was 
again used to compute o(r). The resulting tr{r) curves match 
well onto the curves at smaller radii and there is little depen- 
dence of f, on p within this radial range. 

Also shown in Figure [3] is the influence radius ring of the 
SMBH, defined as the root of the equation 
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Fig. 3. — Relaxation time vs. radius at tlie Galactic cent er. Curves on 
the left (green filled region) assume the mass density of eq. jlOt with three 
normalizations, such that the total (distributed) mass within 1 pc is M» = 
M^,(< lpc)/M,;, = (0.5,1.0,1.5) X 10'*; the = l.Ox lO*" model is shown 
as the heavy line. Curves on the right (blue filled region) assume p « f"''^ at 
all radii with various values for the normalizing density at one parsec, po = 
(0.75. 1.5,3) X l()^M,.)pc"^ the po = 1.5 x lO^Mgpc"' (preferred) model is 
shown with the heavy line. Vertical tick mark is the SMBH influence radius 
computed using c(r) from the preferred model. 

Using the preferred density model with p = 1.5, one finds 
rinfl ~ 2.5 pc. The relaxation time at r w rinfl is a reasonable 
estimate of the time scale over which gravitational encounters 
can change the gross properties of the core. Figure |3] shows 
that the relaxation time at 

The estimates of tr in Figure |3] assume a single population 
of IMp; stars. In reality the Galactic center contains a range 
of mass groups, each of which might have a different veloc- 
ity distribution. Given the relatively long relaxation time, a 
natural case to consider is a distribution of stellar masses, all 
of which have the same velocity dispersion at every radius. 
Equation ^ becomes in this case 
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and N{m)dm is t he number of s tars with masses in the range 
m to m + dm (e.g. ' MerrittI 12004 '). Since equation ((S]) was de- 
rived from the diffusion coefficient ((Av||)") describing grav- 
itational scattering, equation ( I15ab is properly interpreted as 
the time for a test star's velocity to be randomized by encoun- 
ters with more massive objects. Making standard assumptions 
about the initial mass function gives m « IM© (Merritt 2004); 
if the density is dominated locally by stellar black holes then 
m mBH ~ lOM©. It has been argued that m may be even 
larger outside the central few parsecs due to giant molecular 
clouds (iPerets et al.ll2007l) . 

3. MAKING A CORE 



The classical model for a core (e. g. iTremaindl 1 997h is a re- 
gion of constant gravitational potential (|) and constant phase- 
space density / near the center of a galaxy. If the central 
potential is dominated by a SMBH, (|)(r) w — GM,/r and a 
constant / translates into a steeply-rising p: 
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inside ^ rinfl dPeeblesll 19721) . inconsistent with the observed 
distribution. 

The Galactic center has a core of size ^ 0.5 pc, smaller than 
rinfl ~ 2.5 pc (Figs. [U [3]) so equation ( fT6l ) would apply in this 
region. Making a core similar to the observed core therefore 
requires a reduction in the value of / on orbits that pass inside 
- 0.5 pc. 

There are of course many ways to do this. Given that the 
number counts do not strongly constrain the form of n(r) 
within the core (Fig. |2|, we choose not to solve the inverse 
problem n — > /. Instead, we focus on two simple core mod- 
els, both motivated by physical arguments, that are consistent 
with the number count data. 

The starting point is a power-law density of stars at all radii 
around the SMBH, p,- oc r^y^ Q < r <°o. For definiteness, we 
normalize the density of this initial (i.e. coreless) model to be 
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and we assume n,(r) p,(r). (The latter assumption will 
be relaxed in §5.) Equation ( fTTI ) is just equation (fT2l i with 
p = 1.5. The normalizing constant is similar to what vari- 
ous authors have derived in the past for the mass densi ty at 1 
pc; e.g. p « 1.8 (Genzel et al. 2003), p « 2.0 (Schode letaLl 
l2007h etc. The exact normalization is not critical in what fol- 
lows: it serves mostly to fix the relaxation time, and as Fig- 
ure [3] shows, t,- in the region of interest is weakly dependent 
on the density normalization. We show below that our "pre- 
ferred" value of p gives a good fit to the observed stellar ve- 
locities at r w rinfl. 

The stellar mass implied by this model inside one parsec 
is ^ 1.6 X lO^Ma,. This is the mass before the core has 
been carved out; for consistency, this mass should exceed 
the dynamically inferred (distributed) mass inside one par- 
sec, - 1 ± 0.5 X IO^Mq (Scho del etai] |2009). and it does. At 
the same time, the proper motion data are consistent with the 
mass implied by the unmodified power law model, so we can 
not robustly infer the presence of a core in the Galactic center 
mass distribution from the proper motion data alone. 

The gravitational potential generated by p, (r), including the 
contribution from the SMBH, can be written 



GM, 



1-^ 
ro \ r 

1 GMq 
2-Y ro 



1- - 



2-7 



(18) 
(19) 



where ro is a fiducial radius, taken in what follows to be 1 
pc. Mo = IO^M^^Mq is the distributed mass inside ro, and y = 
1.8. The isotropic distribution function fi{E) corresponding 
to the pair of functions (p, , (|)) can be derived numerically from 
Eddington's (1916) formula: 



ME) 



1 1 



dpi d^ 
m VSti:^ dE Ji^(r) # - 



(20) 
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where E — + (|)(r) is the energy per unit mass of a 
star. (Note that the subscript / refers here to "initial," not 
"isotropic") Isotropy is a reasonable assumption, although 
as we will show, some reasonable models for a core imply 
substantial anisotropy. 

3.1. Core origins 

Before proceeding, we consider possible mechanisms for 
the formation of a parsec-scale core at the center of a galaxy 
like the Milky Way. 

1. A binary supermassive black hole. In giant el- 
liptical galaxies, cores are often attribu ted to election of 
stars by a pre -existing binary SMBH (Faber et aP 119971 : 
iMilosavlievic & M erritt 2001), and possibl y to gravitational- 
wave recoil after binary coalesce nce (iBovlan-Kolchin et"an 
120041; iGualandris & Merritt! |2008|) . This model naturally ex- 
plains the sizes of many observed cores - comparable to the 
influence radius of the (single) observed SMBH - if it is 
assumed that the galaxy grew through at least one "rnajor 
merger," with comparably-massive SMBHs (Merritt 2006ah . 
However the binary SMBH model seems less relevant to a 
disk-dominated system like the Milky Way, which may never 
have experienced a major merger Based on a standard ACDM 
cosmological model, the probability that the Milky Way has 
avoided accreting a galaxy with halo mass 1/4 that of the 
Milky Way or grea ter, since a redshift of z = 2, is ^ 30% 
jMerritt et al.ll2002h . The most recent major merger is likely 
to have occurred 10-12 Gyr ago, around the time of formation 
of the thick disk (Wyse 2001). 

2. Inspiral of multiple, smaller black holes. A single 
intermediate-mass black hole (IMBH) of mass ^ IQ^Mq, spi- 
ralling in against a pre-existing stellar density cusp, cre ates a 
core of radius 0.05 - 0.1 pc (iBaumgardt et al.ll2006l) . Re- 
peated inspiral events could create a larger core, although the 
displaced mass would increase at a less than linear rate with 
the number of inspirals. Nevertheless, some models postulate 
one such event every ^ 10^ yr dPortegies Zwart et al.ll2006l) : 
if so, more than one IMBH would probably be present at any 
given time in the inner parsec. 

3. An enlarged loss cone. Gravitational encounters drive a 
mass flux of ~ M,/fr(rinfl) into Sgr A*. The core that results 
from this diffusive loss process is very small: its size is com- 
parable to the radius of the capture sphere - either the tidal 
disruption radius, « 10^^ pc, or the Schwarzschild radius, 
'"Sch ~ 10^^ pc. The core is small because the depleted orbits 
are continuously resupplied by diffusion from orbits of larger 
angular momentum and energy. If there were some way to 
transfer a mass in stars of ^ M, into the SMBH on a time 
scale ^ f, - say, a crossing time - the resulting core would 
be much larger This could happen if the NSC were apprecia- 
bly triaxial, even if only transiently, since many orbits near a 
SMBH in a triaxial cluster are "centrophilic," passing arbitrar- 
ily close to the SMBH after a finite tim e (some multiple of the 
crossing time) jPoon & Merrittil2001h . The size of the resul- 
tant core is determined by a number of factors, including the 
degree of non-axisymmetry and the populati on of the various 
orbit famihes, but it could be of order ^ rinn jMerritt & PoonI 
l2004h . 

4. Localized star formation. Stars might form only, or pref- 
erentially, beyond a certain radius from SgrA*, resulting in a 
low density inside this radius. This possibility is discussed in 
more detail in ^ 1 . 

5. Feedback in active nuclei has also been proposed as a 



rapid core-formation mechanism (iPeirani et al.ll2008h . 

3.2. Isotropic core 

A simple way to lower the density of stars near the SMBH 
is to set 

0, E<Eb. (21) 

Such a truncation leaves the velocity distribution isotropic. 
The cutoff energy can be expressed in terms of a cutoff radius 
rh where 

Eb~^{rb)- (22) 

The configuration-space density after removal of the most- 
bound stars is 

p(r)«p,(r), r>rh (23a) 

k4V2% [ dEf{E)y/(^{r)-E, r<rb. (23b) 

The "approximately equal" sign in these expressions is due to 
the fact that removal of the most-bound stars causes a change 
in the gravitational potential. We ignore that complication in 
what follows, i.e. we replace « by = in equations ( |22] |. ( |23] l 
and similar expressions below. The resultant error in p is at 
most a few per cent, since (|) is dominated by the SMBH at 

^ f'mn, and by the stellar potential due for r ^nfl- 

At small radii, r rb ^ ri^n, the density is 

p^4V2Ky/-iS}{r) [ dEf{E)ocr-^/^. (24) 

In spite of the zero-density hole in phase space, the 
configuration-space density diverges, mildly, toward the 
SMBH. This demonstrates that an isotropic / is not consis- 
tent with a strictly flat core, much less with a central dip in 

P- ^ 

Figure|4]shows space and projected density profiles for var- 
ious values of r/,. The observed number counts are reasonably 
well fit by a model with 

rb w 0.5 pc. (25) 

This value of is close to the core radius derived above from 
the number counts. In other words, in these simple core mod- 
els, rb ~ r core- 

The strange character of these (idealized) models - a zero 
phase space density and a non-zero configuration space den- 
sity - implies some other strange properties. In the Ap- 
pendix, the distribution of orbital elements for stars passing 
near the SMBH is derived. In spite of the isotropic /, the or- 
bits near the center are very eccentric, with an expected size 
a « Q.2rb ~ 0.1 pc. 

3.3. Anisotropic core 

Figure [T| suggests that n{r) may be decreasing toward 
SgrA*. As our second model, we create a core by removing 
all stars on orbits that pass within the sphere r — rb around 
the SMBH. Such a model is a crude description of what hap- 
pens when a binary SMBH ejects stars: in this case, rb would 
be approximately the binary semi-major axis. Collisional de- 
struction of (all) stars that pass within a distance rb of the 
SMBH would also result in such a truncation of /. 

Since we are preferentially removing stars on low-angular- 
momentum orbits, the velocity distribution in the core will be 
anisotropic, biased toward circular motions. 
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r (pc) R (arcsec) 

Fig. 4. — "Isotropic" core models, created by setting / = at low (bound) energies, £<£;, = starting from a power-law model in the density, p « r"'^ 

(shown as the dashed line in both panels). Left panel shows space densities for ri, = (0.1,0.2,0.5, 1,2) pc. Right panel shows surface densities of the same 
models, compared to the number-count data from Figure[T] The vertical normalization of the models in this panel was chosen arbitrarily . 
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Fig. 5. — "Ansotropic" core models, created by setting / = on orbits that pass below ri,. Left panel shows space densities for r;, = (0.02,0.05,0.1,0.2,0.5) 
pc. Right panel shows surface densities of the same models, compared to the number-count data from Figure[T] Blue (thick) curve is marked for comparison with 
the same model in Figure[6] The vertical normahzation of the models in the right panel was chosen arbitrarily. 



We define Eh = ^{rh) as before, and Ec is the energy (ki- 
netic plus potential) of a test particle on a circular orbit of ra- 
dius r/,. If the core is small, ry <^ ri^n, then Ec « E[,/2. Start- 
ing again from a power-law density profile, equation ST% . and 
an isotropic distribution function, equation (|20|) . we now set 
f — when either 



:2rliE- 



-Eb); 



(26a) 



(26b) 



Jb{E) is the specific angular momentum of an orbit with en- 
ergy E and periapse ry. The configuration-space density is 
now zero inside r^: 



p(r)=0, r<rfo. 



(27) 



Outside ry, the density is lowered at all radii due to the ab- 
sence of low-7 stars: 



271 

p '" = — 



Eo(r) 



dEf,{E) 



v/2[£-^(r)]-y2/r 



= 4%/27i:( 1 



1/2 



Eo(r) 



dEf,{E)^E-Eo. 



Here 



(28) 



(29) 



is the minimum energy of orbits that pass through r without 
also passing below r^. 
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r (pc) R (arcsec) 

Fig. 6. — Kinematics of the models of Figure|5] Left panel: anisotropy profiles. Right panel: projected, radial (thick lines) and tangential (thin lines) velocity 
dispersions in the plane of the sky, compared with the observed Or (filled circles) and <5j (open circles) from Schoedel et al. (2009). 



Figure |5] shows density profiles for anisotropic core models 
with various ri,. The observed number counts are reasonably 
well fit by a model with 

rfowO.lpc. (30) 

Note that a given value of r/, produces a larger core than in the 
isotropic core models, since a larger region of phase space is 
affected. 

Defining o, and 0( as the Id velocity dispersions in the ra- 
dial and tangential directions, one finds 



pa? = ^ (l - '^f' l^^dEME) {E-Eof\?U) 




for r > ri,; at smaller radii p = 0. The anisotropy, defined in 
the usual way as p = 1 — af/a^, is then 

3 rl JgdEME)iE-E,)VE^o 

~^^'-d JgdEME){E-Eo?^' 

and is manifestly negative, i.e. o, > O, . 

Figure|6]plots anisotropy profiles for the models in Figure|5] 
as well as projected velocity dispersion profiles in the plane 
of the sky, compared with the observed, proper-motion-based 
velocity dispersions from Schodel et al. (2009). The same 
values of rj, that give a good fit to the number count data, 
also appear to fit the proper motion data reasonably well. The 
tangential anisotropics predicted by these value of - which 
peak between ^ 2 and ^ 5 arc sec - are consistent with what is 
observed, although the statistical error bars are so large that no 
clear discrimination between theoretical models can be made. 

This figure also shows that the adopted density normaliza- 
tion, equation (TU . yields velocities that are consistent with 
the proper motion data. 

4. EVOLUTIONARY MODELS 



The distributions of stars illustrated in Figure [T] and in the 
simple phase-space models of Figures 2] and ID are differ- 
ent from the collisionally relaxed distributi ons normally as- 
sociat ed with stars around a SMBH dBahcall & Wolj|1976l 
IT977h . This may be a consequence of the relatively long re- 
laxation time at the Galactic center (Figure[3]). A key question 
is whether models with a core can survive on Gyr time scales. 
In order to investigate this question, in this section we con- 
sider time-dependent models for the phase-space distribution, 
starting from initial conditions like those discussed in ^ 

Models that are initially isotropic will remain close to 
isotropy as they evolve. Given such initial con ditions, 
the is otropic, orbit-averaged Fokker-Planck equation (iHenonI 
119611) should provide a good description of the evolution. 

Initially anisotropic models will evolve both in /- and E- 
space. Here, we make the approximation that the evolution 
can be divided into two sequential phases: "fast" evolution in 
J, on a time scale <C tr, followed by "slow" evolution in E, on 
a timescale ~ f^. The justification follows from an argument 
made originally by Frank & Rees (1976): Diffusion of stars 
into a point-mass "sink" is dominated by scattering of low-7 
orbits, on a time scale 

fe « e(r)2f,(r) (33) 

where 6 « (r^/ r) '/^ is the angle within which a star's velocity 
vector must lie in order for its orbit to intersect the central 
sink - in our case, the edge of the low-density core at r « rj,. 
The square-root dependence of 9 on fe reflects the fact that 
evolution in 7 is a diffusive process. Thus 

te « -tr (34) 
r 

implying the separation of time scales for stars at r ^ rj,. 

In addition to being physically motivated, this approxima- 
tion allows us to treat the evolution in J and E via differential 
equations with just one space dimension, and to formally sep- 
arate the "isotropization time" from the time scale for changes 
in the radial (i.e. energy) distribution. 

The separation of time scales breaks down for certain orbits, 
e.g. the lowest-energy orbits with E <^ Eh, and a full, /(£,/,?) 
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Fig . 7 . — Evolution of the density of the "anisotropic core" model with = 0. 1 pc due to diffusion in J. Curves show space (left) and projected (right) densities 
at times (0,0.1,0.2,0.5, 1,2,5) x lO' yr. Line thickness increases with time. Other symbols are as in Figure|5] 
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Fig. 8. — Evolution of the space (left) and projected (right) velocity dispersions for the same "anisotropic core" model illustrated in Fig. [7] Solid (dashed) 
curves are radial (tangential) velocity dispersions. Line thickness increases with time; t = (0,0.1,0.5,5) x lO' yr. Symbols have same meaning as in Fig. [6] 



treatment could certainly be justified. Here we note only that 
our initial conditions are somewhat arbitrary, and that changes 
in the /s-dependence of / that would otherwise occur during 
the "fast" evolutionary phase could be seen as establishing a 
slightly different f{E) at the start of the "slow" evolutionary 
phase. 

As the phase space density evolves, the configuration-space 
density p (r) also changes, as well as the velocity dispersions 
a,-(r), Ot{r). Our main constraint on these models is that the 
distribution of stars, after ^10 Gyr, be consistent with the 
observed distribution of late type stars at the Galactic center. 
It is also reasonable to require that a model which matches 
the data now, not be in such a rapid state of evolution that it 
would quickly (in a time <C 10 Gyr) evolve to a very different 
form. 



4.L Evolution in J 
4.1.1. Equations 

The orbit-averaged Fokker-Planck equati on describin g 
changes in / due to diffusion in /-space is (e.g. lCohiilllQTQh 



dN _ d 
dt dll 



(35) 



Here 1{, = J^/J^ is a scaled angular momentum variable, 
Jc{E) is the angular momentum of a circular orbit of energy 
E, N{E,!1{,) = 4%^P{E,!li,)J^{E)f{E,ll) is the number den- 
sity of stars m{E,^) space, P{E, ^ ) is the radial period of an 
orbit, and {D^,D^^^ are the angular momentum diffusion 
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coefficients: 



Z)^(£,^) = -167t2^r2(£) f^fl-^]Fi{E,r), (36a) 



1671^ f dr \ 



Fo{E) 



2 2 



+ 3^FiiE,r) 



F2{E,r) L 



with 



FoiE)=4nT dE'f{E'), 



fi(£,r)=47ir / dE'f{E' 



F2(£,r)=47ir / dE'f{E' 



E 




(36b) 



(37a) 



(37b) 



(37c) 



and r = 4%G^m^\nA. In the expressions (|36] |. the inte- 
gration interval is the radial range from periapsis to apoap- 
sis. Definitions for subsidiary variables can be found in Cohn 
(1979) whose notation is adopted here. Following Shapiro & 
Marchant (1978) and Cohn & Kulsrud (1978), the angular- 
momentum-averaged phase-space density / that appears in 
equations ( l37b is defined as 



m 



d%.f{E,l{.] 



(38) 



Because we are ignoring changes in E, the function f{E) 
does not change with time, nor do the diffusion coefficients 

{Dj^,Dj^j^}. _ 

The practice of some authors (e.g. lMilosavlievic & MerrittI 
|2003) of approximating the 2^ -diffusion coefficients by their 
limiting values as ^ ^> is not followed here. 

Equation (|35] | was advanced in time numerically using the 
NAG routine d03pcf. Initial conditions were f{E,J) corre- 
sponding to the anisotropic core models described in the pre- 
vious section, with various values of r/, (e.g. Figs.|5] |6ll. 

Scaling of the Fokker-Planck models to physical units of 
length and mass is fixed by the power-law density model, 
equation ( fTTl ), used to generate the initial conditions. Since 
the Fokker-Planck equations are orbit-averaged, the relevant 
time unit is the relaxation time. In what follows, times will 
be expressed in years, based on a scaling that assumes a re- 
laxation time given by equation dlSab with m = IMp, and 
InA = 15. If m and InA have different values than these, the 
times given below should be multiplied by 



- /InAX- 



IMqJ \15 J 



4.1.2. Results 



(39) 



Figure |7] shows the evolution of the space and projected 
densities for the anisotropic core model with r/, = 0. 1 pc (blue 
curves in Figs.|5]|6]l. The central configuration-space "hole" is 
rapidly filled as the low-angular-momentum orbits are repop- 
ulated; by a time of ~ 10^ yr, the central density is essentially 
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Fig. 9. — Evolution of the projected anisotropy parameter (s) defined in 
the text, computed over the projected annulus 1" < S < 10". Curves are from 
integrations starting from the five initial models illustrated in Figs.|5]and|6l 
increasing line thickness denotes increasing values off-;,, the initial truncation 
radius, from 0.02 pc to 0.5 pc. Dashed line is the anisotropy at the Galactic 
center as computed from the proper motion data of Schodel et al. (2009), 
with the 90% confidence interval shown as the hatched (blue) region. 

unchanging inside ^ r\j. Because (by assumption) there is no 
diffusion in E, the central hole in phase- space remains, and 
so p(r) evolves asymptoticall to the ~ r^'/^ form demanded 
by an isotropic /(£) with a low-energy truncation (Figure|4|. 
The surface density evolves much less, since even at the (pro- 
jected) center, the surface density is dominated by stars in the 
power-law envelope. Thus, this model is a reasonable fit to 
the observed number counts either at f = or at later times: 
the refilling of the low-/ orbits does not greatly affect the ob- 
served densities. 

Figure [8] illustrates the evolution toward isotropy in the 
same time integration. Initially there is a strong velocity 
anisotropy near the center due to the lack of eccentric orbits 
(Figure|6|. However the refilling of the low-7 orbits increases 

at r r/, on a time scale of ^ 1 Gyr and the core is essen- 
tially isotropic thereafter. 

Proper-motion data from the (projected) inner parsec 
of the Milky Way indicate a slight degree of anisotropy 
(Scho del et al.ll20()9h . We define an averaged anisotropy pa- 
rameter («) as 



(«)^1-. 



(40) 



where a/; and are the radial and tangential velocity dis- 
persions in the plane of the sky and the () denote number- 
weighted averages over some radial range. Adopting 1" < 
R< 10" for this range, the late type stars near the Galactic 
center have 

(«) -0.124°f(f5 (41) 



where the (90%) confidence intervals were derived via the 
bootstrap. Figure |9] compares the observed value of (®) with 
the values predicted by the evolving models. The 90% ob- 
served upper bound, (s) w —0.026, is almost consistent with 
isotropy, which is the asymptotic state of the time-dependent 
models; while the observed lower bound, («) w —0.23, is al- 
most as low as the initial anisotropy of the most extreme core 
model considered here. Thus, the Fokker-Planck models re- 
main consistent with the observed degree of anisotropy over 
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essentially the entire time interval and for a wide range of ini- 
tial conditions. 

Figure |9] also gives an indication of how the time to estab- 
lish isotropy in an initially anisotropic core varies with the 
size of the core. Larger values of ri, imply both a higher ini- 
tial anisotropy, and a longer time scale for the establishment 
of isotropy. 

4.1.3. Summary 

Cores formed by the exclusion of small-periapse orbits 
evolve toward isotropy on a ^ Gyr timescale. This evolution 
does not produce great changes in the observable properties 
of the core, either in the density or the velocity dispersions. 
The anisotropy observed at the Galactic center is consistent 
with the evolving models at both early and late times. 

4.2. Evolution in E 

Evolution in E-space drives / toward the quasi-steady-state 
form 

p^r-y' (42) 

at £■ ^ — GM./rinfl and r <, nnfi, o n a time scale that is 
roughly the relaxation tim e at ^nfl ( iBahcall & Wolj 119761 : 
iLiphtman & Shapirolll977h . Equation ( l42b corresponds to a 
zero net flux in E-space near the hole. In reality, loss of stars 
into the SMBH implies a non-zero flux, causing a gradual 
evolution (expansion) of the cluster, although without much 
change in the form of p(r) (e.g. S hapiro & M archanti ll978l: 
IMurphv etal.l[T99Tl) . We ignore that complication here since 
the density in our models near the SMBH remains far below 
that of the quasi-steady-state models at most times, implying 
a very small flux into the SMBH. 

Many authors have explored quasi-steady-state solutions to 
this equation and to the more general equations that allow for 
a dependence of / o n orbital angula r momentum and stellar 
mass (as reviewed b v lMerrit3l2006bl) . After a. finite time <^ ty, 
the form of f{E) will still reflect the initial conditions. Quin- 
lan (1996) emphasized this in the case of stellar systems with- 
out central SMBHs. Freitag et al. (2006) explored in a limited 
way how the structure of Galactic center models depends on 
the assumed initial density profile. They considered initial 
profiles p ^ r^'' with y as small as 0.75 (steeper, i.e. closer 
to the asymptotic Bahcall-Wolf form, than the steepest initial 
profiles considered here). Freitag et al. found that a time of 
order fr('"infl) is required to erase details of the initial condi- 
tions. Based on Figure |3] that time is ~ 20 Gyr. 

4.2.1. Equations 
The orbit-averaged, isotropic Fokker-Planck equation is 
dN d f df \ 

Dee{E) ^ 16%^T X (43b) 
q{E) j%E'f{E') + j^dE'q{E')f{E') , 

De{E) = 167i:^r j dE'p{E')f{E') (43c) 

(e.g. Cohn 1980; Spitzer 1987). Here N{E) = 4K^p{E) 
is the number density of stars in energy space and p{E) ^ 

4%/2 / ^E - (j)(r)t/r = dq/dE is a phases-space volume el- 
ement. The functions /, De and Dee are understood to de- 



pend on time; as above, the gravitional potential is assumed 
to remain fixed, as do the functions p and q. 

Equations (l43T l were solved using the NAG routine d03pcf . 
The E-space flux was set to zero at the inner boundary of the 
energy grid, as justified above. Two sorts of initial conditions 
were considered: (1) f{E) corresponding to an "isotropic 
core" model; (2) ^-averaged /'s from the final time steps 
of the 7-integrations described in the previous section, which 
started from "anisotropic core" models. (The final /'s from 
these integrations were almost precisely isotropic; cf. Fig- 
ure |9]) In both cases, the initial conditions will be labelled in 
terms of r^. 

4.2.2. Results 

Figure[TO]shows the evolution of the "isotropic core" model 
with rt = I pc. Diffusion in /s-space causes stars to gradually 
occupy orbits of lower (more bound) energies. However even 
after 10 Gyr - roughly the relaxation time at rt, (Figure |3]l - 
/ and p are still far from their steady-state forms at low en- 
ergies/small radii. The Bahcall-Wolf solution is only reached 
after a time that is roughly twice as long. 

The lower panel of Figure[TO]highlights an interesting coin- 
cidence. All of the models considered here have (by assump- 
tion) a density that obeys 

p-r-'-* (44) 

at large radii - the observed dependence of the density of old 
stars on radius beyond ^ 1 pc. But this is essentially the same 
slope as in the steady-state Bahcall-Wolf profile, p ^ ^-i-^s 
which is the asymptotic form of p(r) at small radii. As long 
as the initial core radius is smaller than ^ ri^n, it follows that 
p(r) will evolve in an approximately self-similar way: the 
core will shrink, while outside the core, p (r) will continue to 
obey p ^ r^'^. Reproducing the observed density profile is 
therefore simply a matter of choosing the appropriate, initial 
value of r/,. 

Which values of ri, give cores of the right size now? Two 
estimates were derived in §2 for the core radius of the old 
stellar population, 0.59 and 0.42 pc. These values are plotted 
in Figure[TT] which also shows core radii as a function of time 
in the evolving models (computed in the same way, i.e. by 
finding the projected radius at which the surface density falls 
to 1/2 of its value at 0.04 pc). The figure suggests that the 
currently observed core is consistent with initial cores having 
sizes in the range 

Ipc <; rcore ^ l-5pc, (45) 

or 2 — 3 times the current value. These values are comparable 
with the SMBH influence radius (Fig.O. 

5. SEGREGATION OF THE MASSIVE REMNANTS 

Old stellar populations contain remnants: white dwarves 
(WDs), neutron stars (NSs), and stellar-mass black holes 
(BHs), the end-products of stars with initial masses 1 — 8M0 
(WDs), 8 - 3OM0 (NSs), and - 30 - IOOMq (BHs). Standard 
assumptions about the initial mass function imply that ^ 1 % 
of the total ma ss of a n old population should be in the form 
of stellar BHs (lAlexan der 2005), although observational con- 
straints on t he BH nuniber de nsity near the Galactic center are 
weak (e.g. Muno et al.ll2005h . 

The stellar BHs have significantly higher masses (~ x 10) 
than either the main-sequence stars or the other types of 
remnant that collectively dominate the total mass density. 
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Fig. 10. — Evolution of the phase-space density (top), configuration-space 
density (middle), and surface density (bottom) of the "isotropic core" model 
with rj, = 1 pc (Fig. |4). Increasing line thickness denotes increasing time, 
t = (0,0.2,0.5, 1,2) X lO'" yr. Dashed lines are the asymptotic forms for / 
and p, i.e. / ~ |£ | '^'^j p ~ r^'^l^ . E, f and Z are in arbitrary units. 
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Fig. 1 1. — Evolution of the core radius in various models. Thick curves: 
"isotropic-core" models. Thin cun'es: "anisotropic-core" models. The 
hatched (blue) region encloses the two estimates of the Milky Way core ra- 
dius, as discussed in the text. 



The BHs should therefore lose orbital energy due to dynam- 
ical friction and congregate around the SMBH. Assuming 
that the total density obeys an expression similar to equa- 
tion ( fTTI i. i.e. p ^ r^^, the time for a IOMq BH on a cir- 
cular orbit to spiral all the way in to the cent er is less than 
10 Gyr for a starting radius inside 4 — 5 pc ("M orrisl 119931 : 
iMiralda-Es cude & Gould 2000). Depending on the assumed 
mass fraction in BHs, and on their initial density profile, the 
mass density of BHs after 5 — 10 Gyr is predicted to match 
or exceed that of the other populations inside ^ 10^^ pc. 
If this occurs, the BHs will undergo gravitational scattering 
from themselves and from the other populations, leading to 
a quasi-steady-state n ^ r^^ density profile in the innermost 
regions, and to a slightly shallower profile in the li ghter mass 
components (Freitag et al. 2006; Ho pman & Alexa nder 2006|; 
[Alexander & Hopmanii200 9). 

All of the studies cited above assumed or derived a total 
mass density that increases as p ~ r^^, 1.3 ^ y 2.3 down 
to ~ 10^^ pc from the SMBH. If instead the dominant pop- 
ulation has a core, the dynamical friction force will increase 
more slowly toward the center inside ^ rcore, implying some- 
what longer inspiral times. 

But there is potentially an even more dramatic way in which 
a core can affect the rate of orbital decay, as we now show. 
The instantaneous frictional force felt by a (massive) test body 
of mass ;7iBH and velocity v is 



a = 



A-KG^mBn9{r)F{v)\\\K 



(46) 



where F{v) is the fraction of stars locally that are moving 
more slowly than v . If the latter are described by an isotropic 
/, then 



p(r)F(v) =4%/27i / dEf(E)JE-^r). (47) 



If in addition / is truncated at energies below Et, - an 
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■'isotropic core" - then F falls to zero for orbits with energies 



(48) 



since there are no stars locally that move slower than v for 
these energies - even if the configuration-space density is 
nonzero. Assuming a circular orbit for the test body, and that 
the orbit lies inside the influence radius of the SMBH, this 
condition becomes 

r^'j- (49) 

with Eh — ^{rb)- Thus, inside ^1/2 the core radius, the fric- 
tional force drops precisely to zero. 

This result is not simply an artifact of the brute-force trun- 
cation of /. Consider a core in which p r^^. In a point-mass 
potential, this density is reproduced by 



/(£)=/o|£r^/2 
and F for a circular orbit is easily shown to be 

2 r(y+i) n 



F{y) = 



V7ir(Y-i/2)7,/2 



Jl/2 



(50) 



(51) 



This function varies smoothly from F = 0.5 at y = 2, the sin- 
gular isothermal sphere, to F = at y = 0.5. 

As an even more general illustration of this effect, we com- 
puted the evolution of a circular orbit in a background density 
described by the broken power-law model of equation ([T]i with 
a = 2. The outer slope was fixed at y = 1.8 and the normal- 
ization was chosen to reproduce the density of the "fiducial" 
model, equation ST7\ . outside the core. Setting the inner slope 
y, to 1.8 gives a model similar to those assumed in the stud- 
ies cited above; while small values of y- correspond to a core. 
The radius r of a test body's orbit decays as 



1 dr 
r dt 



'dJ 
|a| I — 
dr 



(52) 



The function F (v) was computed using the expression 
(iSzell et al.ll2005l) 



Fiv)^l~- fd^'%x 
pJE di?' 



1 



tan 



' v/V2 ^ 



(53) 



Figure [12] shows the trajectories of lOM© BHs starting from 
a distance of 4 pc, assuming various values for the inner den- 
sity slope; in each case we set the core radius parameter ro 
to 0.5 pc. As predicted, the rate of orbital decay begins to 
slow when r ro for small values of y, and the decay es- 
sentially stalls, at r < ro/2 pc, when y,- = 0.5. We stress that 
the configuration-space density is nonzero at all radii in these 
models, and in fact increases monotonically toward the cen- 
ter; the precipitous drop in the frictional force is due to the 
lack of low-velocity stars in the core when y,- is small. 

Another consequence of a core is that the dynamical friction 
force along an orbit is not so strongly peaked near periapse, 
which is a necessary condition for an orbit to circularize. In 
the presence of a core, the distribution of BH orbits will there- 
fore remain more nearly isotropic (assuming that it starts out 
isotropic). This, combined with the cutoff in the dynamical 
friction force at r w ro, suggests that the evolved BH density 
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Fig. 12. — Trajectories of lOMy, BHs as they spiral in to the Galactic center 
on circular orbits, starting from a radius of 4 pc. The assumed background 
density is a power-law, p « ' at large radii, with an inner core of radius 
ro, as described by equation (T) with a = 2; ro was set to 0.5 pc (dashed/blue 
Une) and the inner power-law slope of the density was varied, as indicated. 



would rise rapidly toward r w ro, then follow ^ r^^-^ toward 
smaller radii, the density law of an istropic population with an 
inner hole in phase space. If the background density is itself 
evolving, as in the models of the previous section, the core 
radius of the population that produces the dynamical friction 
force will decrease with time, causing the radius of peak BH 
density to also migrate inward on the same time scale. 



5.1. Equations 

As in previous sections, the evolution of the distribution 
of stellar BHs was followed using the orbit-averaged Fokker- 
Planck equation. Define /b//(£', 2^ ) to be the number density 
in phase space of a population of massive objects (BHs), of 
individual mass mBH- We assume that the associated mass 
density Pbh is small compared with the total mass density p 
due to stars and (less massive) remnants. Alexander & Hop- 
man (2009) show that the limiting density ratio for ignoring 
the self-interaction of the BHs is Pbh/P* < ot*/wbh ~ 0.1. 
This condition is violated at late times in some of the integra- 
tions described below; this effectively defines the maximum 
time at which the solutions are valid. 

The orbit-averaged equation describing the evolution of /bh 
as the massive objects experience dynamical friction against 
the background of less massive objects is 



(54) 



where Nm{E,:K.) = 4%^PiE,Ji,)J^{E)fBH{E,%.) is the num- 
ber density in (E,^) space as before. The diffusion coeffi- 
cients depend on the (possibly time-dependent) distribution 
of low-mass objects. Let fi{E,^,t) be the phase-space num- 
ber density of stars with mass m,, m, <^ /«bh- Then (e.g. 
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lTakahashi|[T997l) 



D£(£,^) = -87t2mBH/,^(£)^m, / —Fu{E,r), 

J V,- 



dr 



(55a) 



D^(E,'K.) = -l6%^mzn1lrl(E)Y^mi j ^ ( 1 - ) 
with 



Fu{E,r) ^4%y dE'fi{E',t) 



1/2 



(56) 



and y= 47iG InA; /, is the angular- momentum-averaged / 
as defined above. 

Defining the phase-space mass densities of BHs and the 
other populations (collectively referred to, henceforth, as "the 
stars") respectively as 



gBH = muHfBH , g^Y^mtfi, 
i 

the evolution equation for the BHs can be written 



(57) 




_8_ 

a 

33^ 



gBHJt / - / dE'g{E',t) — ^ 



gBH'K.r, 



.^/^(i-4 

Vr \ 



dE'g{E',t) 



E' -iSf 



1/2- 



This equation was integrated forward using an explicit 
scheme with second-order space derivatives and first-order 
time derivative. Two choices were considered for the stel- 
lar phase-space density g: (1) a time-independent model with 
an "anistropic core"; (2) a time-dependent, isotropic model in 
which g evolves according to equation ( fSSl l. In both cases, the 
phase-space density of the BHs was assumed to be the same 
as that of the stars at f = 0. 

If the mass density in stellar BHs becomes comparable to 
that of the stars at any radius, equations dSST l are no longer 
valid, since the BHs will begin to feel perturbations from each 
other and because the stellar distribution will be affected by 
heating from the BHs. Whether, or when, this occurs depends 
on the assumed initial normalization of ^bh- In plots that fol- 
low, BHs were assumed to be a fraction 10^^ of the total mass 
density initially. 

Results in this section are expressed in years, assuming that 
'Mbh = 1OM0. Times can be scaled to different values of 
men using the simple proportionality of the dynamical fric- 
tion force on wbh (eq.l46]l. 

5.2. Results 

Figure [T3] shows the evolution of the BH mass density 
assuming a fixed stellar background, corresponding to an 
"anisotropic core" model with ri, — 0.2 pc. The initial BH 
distribution is likewise anisotropic. As time progresses, the 
density of BHs drops at a radius of a few parsecs and rises in- 
side ^ 1 pc; BHs accumulate at energies near the core due to 
the falloff in the dynamical friction force there, as discussed 
above. At late times ( ^ 5 Gyr), the BH velocity distribution is 
slightly biased toward circular motions beyond ^ 1 pc, and to- 
ward radial motions inside ^ 1 pc; as noted above, dynamical 



r (pc) 

Fig. 13. — Evolution of the density of stellar BHs assuming a fixed back- 
ground density, shown by the thin (blue) line. Curves show the density at 
times (0,1,2,4,8,16) Gyr. 



friction in the presence of a core does not efficiently circu- 
larize orbits. To a reasonable approximation, the phase-space 
distribution of the BHs at late times is isotropic, and because 
the stellar distribution has not been allowed to evolve, the BH 
(5 5)density remains zero on orbits with energies below E^ ~ Eh/2. 
The result is a p^// ^ ^-0-5 (.^gp at r ^ 2ri « 0.4 pc. Note that 
the density of BHs just approaches that of the stars at the final 
time (16 Gyr) in this integration. 

Allowing the stellar distribution to evolve is more realistic. 
Figures fl4lfT5] show the results of two such integrations. In 
Figure [141 he initial distributions of stars and BHs were gen- 
erated from an "isotropic core" model with rt, ~ 2 pc. The 
stellar distribution was allowed to evolve according to equa- 
tions (l43T l. yielding a time-dependent g{E,t) which was in- 
serted into equations ( fSSl l at each time step to compute the 
diffusion coefficients acting on gsH- In these integrations, the 
stellar core shrinks, on a time scale that is ^ 10 times longer 
than the time scale for the BHs to accumulate around the core. 
As a result, the BHs "follow" the stellar core inward. Their 
density does not rise so steeply toward the core as in the in- 
tegrations with fixed stellar density (Figure [13]) since the ra- 
dius at which they would otherwise accumulate changes with 
time. Figure[T5]shows the mass enclosed within (0. 1 , 0.3, l)pc 
vs. time for both components, in a second integration start- 
ing from ri, ~ I pc. Also shown for comparison is the mass 
in BHs estimated by Miralda-Escude & Gould (2000) to lie 
within the central parsec after 10 Gyr, assuming dynamical 
friction against a fixed stellar background; and the mass in 
BHs estimated by Hopman & Alexander (2006) to lie within 
0.1 pc, based on their steady-state multi-mass Fokker-Planck 
solutions. 

The feature that we wish to emphasize here is the sensitiv- 
ity, in our models, of the final density in BHs to the elapsed 
time. Figures[T4]and[T5]suggest that it would be unjustified to 
assume that the stellar BHs have reached a steady-state den- 
sity by now at any radius inside ~ 1 pc. This is even more true 
if star formation has been an ongoing process in the nuclea r 
star cluster (e.g. ISerabvn & Morrislll99^ ; 'Figere t al.ll2004D . 
since the mean age of stars and remnants may be much less 
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Fig. 14. — Evolution of the number density of stellar BHs (thick/black 
curves) assuming an evolving background (stellar) density (thin/blue curves), 
starting from an "isotropic core" model with r;, = 2 pc in both components. 
Times shown are (0, 1,2,4,8) Gyr. 
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Fig . 15. — Evolution of the enclosed mass in stellar BHs (das h-dotted lines) 
and stars (solid lines) in an integration like that of Figure [14] with rj, = 0.1 
pc. Curves show mass enclosed within (0.1,0.3, 1.0) pc assuming that stellar 
BHs have initially 1% the mass density of stars at each radius. The hatched 
regions show the estimates of Miralda-Escude & Gould (2000) for the mass 
of stellar BHs within 1 pc after 10 Gyr, assuming dynamical friction against a 
fixed stellar density cusp; and of Hopman & Alexander (2006a) for the mass 
of BHs within 0. 1 pc based on steady-state Fokker-Planck solutions. The 
vertical width of each hatched region corresponds to an (arbitrary) factor two 
in mass. 

than 10 Gyr (EB- 

This point is made more forcefully in Figure [16] which 
compares stellar and BH densities at r = 0.01 pc in three 
models with different initial core sizes. The figure shows that 
when the evolving stellar core reaches a size consistent with 
the observed core, the density in stellar BHs can be substan- 
tially less than that of the stars at this radius. In the steady- 
state models, the BH density meets or exceeds that of the stars 
at radii of lO^^ - lO"' pc (e.g. lFreitag e"tar..,2006. Fig. 10), 
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Fig. 16. — Density of stellar BHs (dash-dotted curves) and stars (solid 
curves) at r = 0.01 pc in three models with cores of different initial size, 
rt = (0.5, 1,2) pc. Horizontal axis is core radius of the stars, which decreases 
with time as the core shrinks. Integration times are 10 Gyr. Vertical hatched 
region indicates the radius of the observed core. 



and gravitational wave inspiral of BHs is dominated by BHs at 
these distances (iHils & Bendeii [19951 : iHopman & Alexandeii 
I2005h . The much lower densities found here have poten- 
tially important implications for the predicted rate of inspiral 
events, as discussed in more detail in 36.41 

In models with initial core radii ^ 1 .5 pc (small enough to 
reproduce the currently-observed core), the density of stellar 
BHs after 10 Gyr becomes large enough that self-interactions 
between the BHs would be significant. A high enough den- 
sity of BHs would also tend to accelerate the relaxation of the 
stellar component. While beyond the scope of this paper, tests 
of these predictions could be carried out via two-component, 
f{E,L,t) models or A^-body integrations. 

5.2.1. Summary 

The time scale for inspiral of IOMq BHs to the center of the 
Galaxy is longer in models with a core than in models with a 
cusp, both because of the lower density of stars in the core, but 
also because the dynamical friction force drops essentially to 
zero at energies near the phase-space truncation energy that 
defines the core. Assuming that the BHs and the stars have 
the same distributions initially, the radius of peak density of 
the BHs tends to "follow" the stellar core as the latter shrinks. 
The density of BHs can remain substantially less than that of 
the stars at small radii (r <^ 0.01 —0.1 pc) even after the stellar 
core has shrunk to its observed size of ^ 0.5 pc. 

6. DISCUSSION 

Topics discussed in this section include the dynamical con- 
sequences of ongoing star formation; predicted rates of stellar 
tidal disruptions and gravitational wave inspirals at the Galac- 
tic center; mechanisms for enhanced relaxation; the effect of 
a parsec-scale stellar core on inspiral of intermediate mass 
black holes; the rate of production of hypervelocity stars by 
collisions involving stellar BHs; and the connection between 
the core at the center of our Galaxy and the cores observed in 
other galaxies. 



16 



David Merritt 



6.1. Star formation 

The MiUcy Way nuclear star cluster (NSC) sits at the cen- 
ter of a kiloparsec-scale Galactic bulge or bar which consists 
mainly of old 10 Gyr) evolved stars. On sub-kpc scales, 
the Milky Way shows evidence for stellar populations with a 
range of ages. Serabyn & Morris (1996) argued that the con- 
ditions of the interstellar medium in the "Central Molecular 
Zone" would lead inexorably to inflow of molecular mate- 
rial and to continuous star formation activity in the central 
~ lO'^ pc. Figer et al. (2004) modelled the luminosity func- 
tion of stars in the inner ^ 50 pc and argued that single-burst 
star formation models could be securely ruled out; they in- 
ferred a nearly constant rate of star formation rate over the last 
~ 10 Gyr. These a nd other studies (e.g. Mezgeret al. 1999; 
iPhilipp et aDll999h suggest that the NSC is not a simple in- 
ward extrapolation of the old bulge, but rather consists of an 
intermediate age population that has been undergoing contin- 
uous star formation since the creation of the Galaxy. 

In a general way, continuous star formation strengthens the 
picture presented here of a nuclear cluster that is less than ^ 
one relaxation time old, by reducing the mean age of stars 
from 10 Gyr to 5 Gyr. 

Of more direct interest is the evidence fo r recent star for- 
mation in the inner ^ O.Spcdl aumar dl2008h . the same region 
where the old stars exhibit a low-density core. The total mass 
of the young stars currently observed in this region is proba- 
bly less than IO'^Mq (Figerl l2008h . making them dynamically 
insignificant. On the other hand, the starburst that created this 
population may have been just the most recent instance of an 
ongoing or episodic process. How would such a "source term" 
modify the evolutionary calculations presented above? 

The answer clearly depends on the accumulated stellar 
mass and on the radial dependence of the star formation rate; 
both are highly uncertain. Here we limit ourselves to answer- 
ing a simpler question: How would a population of stars, 
formed initially in a disk, evolve against the background of 
a pre-existing stellar core? We specifically ignore interac- 
tions between the disk stars, and assume that the background 
stars are fixed in their distribution. More detailed calculations, 
in which both populations are allowed to evolve, will be de- 
scribed in a subsequent paper. 

If the surface density of the young stellar disk is 
E(r) o<: r^", the distribution of orbital energies is dN/dE = 
{dN/dr){dr/dE) °^ \E\"^^. To simplify the evolutionary cal- 
culation, we suppose that the orbits of the young stars "ran- 
domize" in orientation and eccentricity on a time scale shorter 
than f,.. (Without trying to justify that assumption quan- 
titatively, we note several mechanisms that might achieve 
this: standard relaxation as in §3; "resonant relaxation" as in 
§6.4; torques from another disk or from a large-scale non- 
axisymmetric bulge component etc.) The young stars can 
then be described by an isotropic / where f{E,t = 0) o<: 
p{E)-^N{E) oc |£|"-V2. The corresponding space density is 
p{r,t — 0) oc r^^^". If the disk is truncated at an inner radius 
Tin, then p(r,f — 0)°^ ^-1/2 inside ^ rin. We assume « = 2 and 
that the disk stars are initially distributed between 0.1 pc and 
0.4 pc. 

For the old stars, we assume a (fixed) /oid(£^) = /o|£^|'* 
which corresponds to a density Poid('") — Po'"^, y — A + 
3/2, r < rinfl. Substituting this expression for f^M into the 
energy-space diffusion coefficients De,Dee of equation ( l43T l 



gives 

DEiE) = -^^foG'MWinA\Ery\ (59a) 

= (1+A)(?"-2A) ^''^'^'"'^"^''^'^''^^^^^^ 

again assuming that the gravitational potential is due to the 
SMBH alone. 
We adopt units for time and energy such that 



4jtr/o[r][£]^ = i 



where F = In A. A natural energy unit is 



[E] = ¥0 = 



GM. 



(60) 



(61) 



which makes the unit of time 

[T]-^ = leTZ^GWlnAfoy]!^. 

In dimensionless variables, equation ( |43] ) for the young stars 
then becomes 



(62) 



dt dE' 



(63a) 



F(E) ^ -E^— (E'^f) (63b) 

^ ' 1+A) 1-2A) dE^ ^ ' 



(1+A)(1-2A) 



where 



^^-(1+A)(1-2A) -1+6A-8A^ (64) 

3-2A ' 2(3-2A) ^ ' 

The steady-state solution is given by setting F{E) —0: 



f{E,t 



(65) 



Choosing A = 1 /4 (y 7/4) implies B = - 1 /4; in other 
words, test particles interacting with stars in a Bahcall-Wolf 
cusp evolve also to the Bahcall-Wolf form. However, other 
choices for A imply different steady-state solutions for the 
young stars. 

Figure [TT] shows time-dependent solutions to equation ( l63T l 
under two assumptions about the background population. In 
both cases, the (fixed) density of the old population was nor- 
malized to give a mass within one parsec of 1.5 X IO'^Mq. 
Figure [TTh assumes Poid °= 

-7/4 (A = 1 /4), a Bahcall-Wolf 
cusp. In this case, the relaxation time at the initial disk radius 
is just a few Gyr, and the young stellar population also reaches 
a distribution close to the Bahcall-Wolf form after 5 Gyr, as 
expected. Gravitational encounters with the old stars tend to 
redistribute the young stars to both smaller and larger radii. 

Figure [T7b shows the case Poid °= r^^/^ (A = —3/4), a 
lower-density core. Now, the relaxation time increases toward 
the center; the distribution of young stars hardly changes in- 
side the original inner disk radius, and the net result of gravi- 
tational encounters is mostly to scatter the disk stars to larger 
radii. The steady-state form, p oc is only gradually 

approached, and only at radii outside the initial disk radius. 
There is essentially no evolution toward a Bahcall-Wolf cusp. 

In both cases, the number of young stars that remain within 
the original outer disk radius after 5 Gyr is a small fraction of 
the initial number: ~ 10^^ in the A = 1/4 case and ^ 0.15 for 
A = -3/4. 
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Fig. 17. — Evolution of f{E,t) and p{rj) describing an isotropic popula- 
tion of "young" stars that formed initially in a disk, and subsequently scatter 
off of the "old" stars, assumed to have a fixed phase-space distribution. Ini- 
tial conditions are shown in bold; subsequent times are (1,2,3,4,5) x lO' yr, 
assuming that the mass in the old population within 1 pc is 1.5 x 10*M,;;. 
Length units were scaled to the Galactic center, (a) /ojj « p^u <x r^'^/'^, 

the Bahcall-Wolf form, (b) « E^^^*, paid °= r^^^*. Dashed lines show the 
steady-state slopes. 

These results, while very preliminary, suggest that ongo- 
ing star formation need not greatly modify the conclusions 
that were arrived at above. As long as star formation occurs 
against the backdrop of a stellar core, the density profile of 
the young stars can remain relatively flat inside the initial disk 
radius. Their mean density is also strongly diluted by encoun- 
ters with the older stars. 

This conclusion may need to be modified in the case of 
the stellar remnants. The ^ 10 — 100 very massive blue gi- 
ant stars in the stellar disks will probably end their lives as an 
equal number of IOMq BHs, with total mass ^ lO^-lO^M,... 
If such star forming events occur once per 10^ yr, then 
IO^Mq — IO^Mq in stellar BHs could accumulate in the cen- 
tral parsec over 10 Gyr Q The expected number of remnants 
will also depend strongly on the form of the initial mass func- 
tion ( IMF); for instance, a "top-heavy" IMF (.Paumard et al.l 
120061) would produce many more stellar BHs per unit of total 
mass. 

6.2. Rates of tidal disruption 

The existence of a core implies a smaller rate of stellar tidal 
disruptions at the Galactic center than in models that assume 
a density cusp. The rate of scattering of stars into the SMBH's 
tidal disruption sphere, r < r,, can be written 

N = J j{E)dE (66) 

where (£) is the number of stars scattered per unit time and 
unit energy into r,. In the models considered here, setting 
the core radius to zero gives p r^^-^ at all radii (eq. [TZb . 

' I thank Tal Alexander for pointing this out. 



TABLE 1 
Stellar tidal 
disruption rates 



rb (pc) 


W(yr-') 





4.6 X lO-'' 


0.1 


2.7 X 10"" 


0.2 


1.7 X lO""* 


0.5 


5.6 X 10-5 


1 


2.7 X 10-5 


2 


2.0 X 10-5 



and this is also approximately the form of the density in the 
single-mass Bahcall-Wolf steady-state solution at r <; rinfl, 
p r^''^^. Both density laws, in turn, are close to a singu- 
lar isothermal sphere, p o= r^^, for which the feeding rate has 
been shown to be 

A. . 4.6 X 10-V-' { " { \ ) " 

V90kms-V \Ax\Q('M,;J 

(67) 

assuming Solar-type stars ( IWang & Merrittll2004t) . 

Carving out a core changes for two reasons. (1) The en- 
ergy integral, equation ( |66] |, now has Eb as a lower limit rather 
than — oo. (2) The diffusion coefficients IX^D.j^.j^ that deter- 
mine the scattering rate at every E (eqs.[36ll are smaller due 
again to the absence of stars with E < Ej,. 

A commonly made approximation (e.g. ' Cohn & KulsrudI 
119781) is to ignore the contribution to F{E) from scattering 
off of stars with energies that are smaller (i.e. more tightly 
bound) than E. In this approximation, the diffusion coeffi- 
cients are not changed, at energies E > Ej,, by the truncation 
of / at Eh, and the only change in comes from the change 
in the lower integration limit for equation ( |66] ). 

Table 1 shows for the Milky Way computed under this ap- 
proximation, for various values of r^, assuming an "isotropic 
core." (Recall that rt, is essentially equal to the core radius 
of the corresponding density profile. Figure H] so that setting 
r/, « 0.5 pc gives a core of roughly the correct size for the 
Milky Way.) A 0.5 pc core implies a tidal flaring rate that is 
almost an order of magnitude smaller than for a coreless cusp. 

This conclusion should be considered extreme since it ig- 
nores changes in / that will necessarily tend to refill the de- 
pleted orbits over relaxation time scales. In principle, the time 
integrations of / in §4 could be used to compute the evolution 
of the tidal disruption rate. We postpone that calculation to a 
later paper, but make a related point here: if / is not in a steady 
state, the standard expressions for the flux into the loss cone 
(like the expr essions that were used to deriv e eq.l67]) are not 
strictl y valid (iMilosavljevic & Merr"it3l2003l: [Merritt & Wan^ 
l2005h . This highlights the need to develop a more complete 
theory of time-dependent loss cones. 

6.3. Enhanced relaxation 

The evolution equations for / and /bh that were solved in 
§4 and §5 were based on a standard, orbit-averaged Fokker- 
Planck treatment of gravitational encounters. The results were 
scaled to physical time units assuming that relaxation is driven 
by perturbers with masses of roughly Mq. 

Both of these assumptions have been questioned, in recent 
papers that argue for more efficient relax a tion near the Galac- 
tic center dHopman & Alexandeii 120061: iPerets et alj l2007h . 
Two important points have been made: (1) The effectiveness 
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of gravitational encounters at inducing changes in orbital an- 
gular momenta is increased for stars orbiting within the grav- 
itational field of a SMBH, r < rinfl, since the orbits are nearly 
Keplerian and they maintain their orientations for many radial 
periods, allowing torques to accumulate linearly with time 
("resonant relaxation"; Rauch & Tremaine 1996). (2) If there 
is a distribution of masses in the scattering population, the ef- 
fective relaxation time is determined by the second moment 
of the mass function (eqs.fTSl). The presence of even a small 
population of very massive objects ("massive perturbers") in 
the region of interest can reduce the effective relaxation time 
considerably. 

We briefly discuss the applicability of these ideas to the evo- 
lutionary models discussed here. 

Resonant relaxation: Resonant relaxation is relevant to the 
timescale for isotropization that was computed in §4. The 
resonant relaxation time is 



1 /M. 

A^(< a) \m* 



PHa) 



(68) 



precess 



jHopman & Alexanderll2006l) where N{<a) is the number of 
perturbing stars, of mass m*, inside the orbit of the test star 
whose semi-major axis is a; P is the orbital period; and fpiecess 
is the time scale above which orbits lose their coherence due 
to precession. If precession is due primarily to the gravita- 
tional force from the stars, distributed spherically inside 
r = a, then 

1 M. , , 

-P{a) (69) 



and 



precess ' 



fRR ~ P{a) 

11, 



A^(< a) nii. 



^2x lO^'yr 



3/2 



(70a) 
(70b) 



independent of the density of perturbers. Hopman and 
Alexander (2006b) note that this time falls below the standard 
relaxation time ([8]) at a distance 0. 1 — 0.5 pc from SgrA* . This 
suggests that resonant relaxation may reduce somewhat the 
time scale for isotropization compared with the values com- 
puted in §4. 

The counter-intuitive result that fRR is independent of the 
density of perturbing stars (eq. I70bb is due to the fact that the 
coherence time fprecess becomes long as becomes small, al- 
lowing even small torques to build up for long times. In the 
core models discussed here, can be essentially zero, and 
other processes would likely begin to d ominate the preces- 
sion rate. For instance, a nuclear bar (e.g. lAlardll2001l) would 
generate a non-axisymmetric component to the gravitational 
potential in its interior, setting an upper limit to the preces- 
sion time and reducing the effectiveness of resonant relaxation 
compared with the expressions given above as — >^ 0. 

Massive perturbers: Equation ( fTSl l says that the rate of 
gravitational scattering by a background population of per- 
turbers scales as 

np{ml) (71) 

where np is the number density of perturbers of individual 
mass nip-, the brackets denote a number-weighted average. 
This expression ignores differences in the velocity distribu- 
tion between the different populations, which is reasonable if 



the relaxation time is long and/or if the perturbers were re- 
cently formed. Perets et al. (2007) noted that the mass spec- 
trum of giant molecular clouds implies that they should dom- 
inate the scattering rate beyond a few parsecs from SgrA*, 
reducing the relaxation time in this region by as much as 
several orders of magnitude. Inside ^ 5 pc, they suggested 
that gas clumps in the circumnuclear gas disk, with masses 
10^ — lO^Mp), might be similarly important. 

Massive perturbers beyond ^ 1 pc would affect the distribu- 
tion o f stars near the SMBH in two distinct ways (IPerets et alj 
|2007|). 1. Deflection of unbound (with respect to the SMBH) 
stars onto radial orbits would fill in some of the phase space 
volume that was evacuated by formation of the core. This 
population would have a spatial distribution n ^ r^^l'^, simi- 
lar to that of the pre-existing core. 2. Three-body interactions 
of field binaries deflected by massive perturbers can create 
a population of bound stars around the SMBH, initially on 
very eccentric orbits. The radial distribution of these bound 
stars will reflect the semi-major axis distribution of the par- 
ent binary population, which is uncertain. The capture rate is 
estimated to be as large as ^ lO^"* yr^' in the inner parsec; 
the accumulated mass could therefore potentially exceed the 
number of stars in the evolutionary models considered here. 
These arguments, while very approximate, suggest that mas- 
sive perturbers could compete with stellar-mass perturbers in 
terms of refilling an evacuated core. 

In a general way, the fact that the Milky Way does contain a 
low-density core implies an upper limit on the effectiveness of 
any relaxation process, particularly those that change orbital 
energies. 

6.4. Extreme-mass-ratio inspirals 

Inspiral of compact remnants (stellar mass BHs, neutron 
stars, white dwarves) into a SMBH is accompanied by the 
emission of gravitational waves with frequencies that will 
be detectable by the Laser Interferometer Space Antenna 
(USA) dSjg urdsson & Rees 1997; iBarack & Cutlej i200l 
lAmaro-Seoane et al.. 2007) . Event rates of these extreme- 
mass-ratio inspirals (EMRIs) are generally computed as- 
suming that the stars and stellar remnants are distributed 
in a relaxed, multi-mass density cus p (Hils & Bender 1995]; 
iHopman & Alexandeiil2005t lHopmanll20 09). In such a cusp, 
the density of solar-mass stars follows n ^ r^^-^ while the 
lOM© BHs have a steeper dependence, n ^ r^^-^^ — r^^. The 
radius at which the mass density of BHs rises above that of the 
less massive objects depends on the choices made for the mass 
function, and (in the time-dependent models) the initial spa- 
tial distributions an d the elapse d time. Ty pically, Prh > P* in- 
side .01 -0.1 pc (iFreitag et al . 2006; H opman & Alexander 
120061) . The EMRI event rate is dominated by BHs (as opposed 
to neutron stars or white dwarves) due to their high masses 
and high mass densities. Most of the signal is contributed by 
BHs inside ^ 10^^ pc (e. g. Hopman & Alexander 2006a, 
Figure 2). 

Here we estimate BH inspiral rates for models of the nu- 
clear star cluster that include a core. A distribution-function 
based approach, similar to what was used above to com- 
pute stellar tidal disruption rates, would require an additional 
Monte Carlo calculation to estimate the probability that a 
star on a loss cone orbit will evade being scattered directly 
into the SMBH before emitting gravitational waves. Instead, 
we follow the more approximate treatments in Hils & Ben- 
der (1995), Hopman & Alexander (2005) and Ivanov (2002) 
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based only on the density profiles of the two components. We 
first derive an expression for the inspiral rate based on the 
observed stellar distribution. Figures 1 and 3. Event rates pre- 
dicted by the evolving, two-component models of §5 are then 
computed. 

Gravitational wave emission is dominated by BHs that are 
scattered into the SMBH from orbits with a < a^nt, where a^nt 
is the orbital radius such that the decay time due to emission of 
gravitational waves equals the time for the BH to be scattered 
in or out of the loss cone by stars or by other BHs. The angular 
momentum of a loss cone orbit is 



Ac 



4GM, 



(72) 



Stars with a <^ rcrft can avoid scattering for a time long enough 
to spiral in via emission of gravitational waves. The inspiral 
rate is approximately 



A^GW 



/BH(a)A^(a) 
fr(fl)ln0 



da 



(73) 



where N{a)da is the number of stars and BHs with semi- 
major axes a to fl + da, /bh < 1 is the fraction of objects which 
are BHs, and = Jc/J\c- 

The gra vitational wave inspiral time for an eccentric orbit, 
J < Jc, is jHopman & Alexanderil2005h 



3 X ^/GM^ M. 
85 niBH 



(74) 



where msH is the mass of a stellar BH. The critical radius is 
defined as the radius where few — he with 



(75) 



the diffusion time into the loss cone. Combining (|74] | and 
dTSl ), and identifying the relaxation time on an orbit with its 
value at r — £3, gives an implicit relation for flcrit- 



85mBH fGM, 



1/2 



fr(flciit) 



,10 



(76) 



the solution to which determines the upper limit to the rate 
integral ( |73T l. 

For the relaxation time, we use the estimate plotted in Fig- 
ure|3]for the Galactic center, which was based on the paramet- 
ric model fit to the number counts. Figure [1] with n oc r^'/^ 
enforced at small radii. This is the flattest central dependence 
consistent with an isotropic / for the stars (§3), and is also 
similar to what is found in the evolving models (§4, 5). The 
relaxation time plotted in Figure|3]can be written 



tr{r) « 1 X 10' V 



m 



-1 



0.2pc 



-1 



r<0.2pc (77) 



where m is defined in equation ( llSbb and accounts for the 
possibility that stellar BHs may contribute significantly to the 
density of scatterers. (Note however that this expression as- 
sumes a particular normalization for the total mass density; 
this assumption will be relaxed below). In our two-component 
model, m is fixed by /bh- 



m _ 1 -/bh + /bhWbh/»^ _ 1 + 99/bh 
m 1 - /bh + /bhotbh /»! 1 + 9/bh 



(78) 




t (Gyr) 

Fig. 18. — Evolution of the BH inspiral rate in the models of Figure ll6l The 
solid part of each curve terminates at the time when the stellar core radius is 
0.5 pc. Dashed (blue) line shows the event rate computed by Hopman & 
Alexander (2()06a) in steady-state, mass-segregated models (cf. their Fig. 2). 

with m the stellar mass; the last relation assumes iiibh = lOm. 
Substituting (iTTl i into ( |76] | gives for the critical radius 



flerit - 0.08 pc ( —] 



-0.4 



OTBH , 



(79) 



For m/mBH = 0.1(0.5)1, Ocrit/pc = 0.2(0.11)0.08. These are 
somewhat larger than the critical radii computed assuming a 
relaxed density cusp (ocrit ~ 10^^ pc) but are still small com- 
pared with the observed core radius of ^ 0.5 pc. 
The event rate ( |73] l becomes 



Afew 



:3x lO^Gyr-'/BH 



'«BH 



-0.4 



jVo.i/10^ 
In 0/5 



(80) 



where A^o.i is the number of stars and BHs inside 0.1 pc; 
the density normalization assumed in making Figure [3j equa- 
tion ([Toll, implies A'o.i « 8.6 x 10^ 

Setting /bh = 0.001 and m « IMp, in this formula gives 
an estimate of the inspiral rate for an unsegregated model 
in which the BHs follow the same density profile as the ob- 
served stars. Not surprisingly, the resultant rate is very low, 
N Ki 1 Gyr^'. If instead /bh is set to the higher values found 
in the steady-state models, i.e. /bh = 0.01(0.05)0.1, equa- 
tion dSO]) gives A^GW ~ 2(9)15 Gyr"'. These rates are still 
1-2 orders of magnitude lower than those in the steady-state 
models, due to the very different total densities assumed at 
the relevant radii; for instance, Hopman & Alexander (2006a) 
find Afew ~ 300 Gyr"' in models that have 0.01 <, /bh ^ 0. 1 
at r « flciit. The magnitude of the difference is at first sight 
surprisingly small given that the steady-state models have to- 
tal densities at 10"^ — 10"^ pc that are orders of magnitude 
higher than here. The reason is that the longer relaxation time 
in the core models implies a ^ lOx larger value of a^nt and 
a ~ lO^x larger volume fro m which inspirals can occur (e.g. 
iHopman & Alexandeiil2005h . 

In the time-dependent models of §5, not only are total den- 
sities smaller than in the steady-state models, but /bh also 
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remains well below its value in a mass-segregated cusp until 
late times. Figure[T8]shows BH inspiral rates computed in the 
same way as above, using the densities of stars and BHs in 
the time-dependent models. In these models, the BHs have 
initially the same spatial distribution as the stars; as the stellar 
core shrinks, the BHs "follow" it inward, increasing both their 
density near the center and their relative density with respect 
to the stars (Figure [TSIl. By the time the stellar core shrinkss 
to its observed value of ~ 0.5 pc, the BH density at r <, cicrit 
is still much lower than in the steady-state models and Afew is 
also correspondingly lower. 

It must be emphasized that these evolutionary models are 
idealized and that one could imagine other, reasonable initial 
conditions that would produce rather different rates after 10 
Gyr For instance, the distribution of BHs might have exhib- 
ited some degree of segregation even at early times. The main 
point to be made here is that in models for the nuclear star 
cluster that contain a core, Nqw could plausibly be as much 
as 10-100 times lower than in models based on a steady-state, 
mass-segregated density cusp. 

6.5. Inspiral of an IMBH 

Intermediate mass black holes (IMBHs), with masses of 
10^ — lO^M©, may form in de nse star clusters through run- 
away mergers of ma ssive stars ( iPortegies Zwart & McMillan! 
l200a iFreitag et al.1 |2006i) . Inspiral of an IMBH into the 
Galactic center is usually modelled assuming that the dynam- 
ical friction force originates in a stellar density cusp with p ^ 
r_2J.5^Jl^_2 (Baumgardt et al. 2006; Matsubayashi et aL 
120081; iLockmann and Baumgardt 2008). In these circum- 
stances, inspiral continues until the IMBH reaches a distance 
from SgrA* such that the enclosed stellar mass is roughly 
equal to ;7Jimbh, or lO^-'pc ^ r <, lO^^pc. 

If there is a pre-exisiting core in the stellar distribution, in- 
spiral would stall at roughly 1 /2 the core radius, or ^ 0.25 
pc in the case of the Galactic center (Figure [T2]l. independent 
of the IMBH mass. Furthermore, as noted above, the orbital 
eccentricity of the IMBH would not be expected to decrease 
strongly during the inspiral. 

Men-itt et al. (2009) and Gualandris & Merritt (2009) noted 
one consequence of a "stalled" IMBH on an eccentric orbit: 
the IMBH very efficiently randomizes the orbits of ambient 
stars, producing a nearly "thermal" distribution of orbital ec- 
centricities, N{< e) ^ e^, on Myr time scales. Merritt et al. 
(2009) postulated stalling radii inside ^0.1 pc in order to ex- 
plain t he observed distribu tion of the S-star orbits (Gh ez et al.l 
120081; iGillessen et"an 1200 9). The somewhat larger stalling 
radii made plausible here suggest that IMBHs might be lurk- 
ing on somewhat wider orbits, roughly the size of the two stel- 
lar di sks at O.lpc ^ r < , 0.5 pc (Paumard et al. 2006; Lu et 
|^09';'Bartko et al. 2009a). As shown by Levin et al. (2005) 
and Berukoff & Hansen (2006), an IMBH at these radii could 
play a role in truncating the stellar disks and scattering disk 
stars onto inclined and eccentric orbits. Those authors as- 
sumed efficient inspiral of the IMBHs, which limited the time 
over which interactions could occur. If orbital decay stalls at 
distances of ^ 0.2 pc, dynamical interactions with disk stars 
could be prolonged indefinitely, potentially resulting in much 
greater changes in the stellar orbits. 

An inspiralling IMBH also ejects stars and stellar rem- 
nants via three-body interactions with the IMBH/SMBH bi- 
nary. Some of these stars receive kicks greater than ~ 



10^ km s^', allowing them to escape into the Galactic 
halo as hyper- velocity stars ( HVSs ) ( Yu & Tremaine 2003). 
Simulations of this process (Levin 2006; Baumgardt et al] 
2006; Lock mann and Baumgardt 2008) also typically assume 
a steep density cusp for the stars. The resulting rate of ejection 
of HVSs increases rapidly with time as the IMBH spirals in, 
peaking when the IMBH reaches its (small) stalling radius of 
^ 0.01 pc, then falling off due to the local depletion of stars. 
The stellar density at this radius is assumed to be > lO^M© 
pc^^ initially, roughly the density implied by the (coreless) 
power-law model of equation (fT2l i. 

In the core models considered here, the rate of production 
of HVSs would be much smaller than in the relaxed models 
due to the lower stellar densities inside ^ 1 pc. In addition, 
the IMBH would stall at a larger radius of ^ 0.2 pc. The den- 
sity at this radius is lO^M© pc^^, resulting in a ~' 10^ times 
smaller rate of escapers than the peak values obtained during 
infall in models with a steep cusp, i.e. 0.1 — 1 Myr^'. On the 
other hand, at these low rates, ejections by the IMBH would 
hardly affect the ambient stellar density and the production 
of HVSs could continue indefinitely at an approximately con- 
stant rate. Production of ~ 10 HVSs would therefore require 
a span of ^ 10—100 Myr; in fact the observed span of travel 
times for escaping HVSs is 2 x 10^ yr (iBrown et al.ll2007l) . 
roughly consistent with this crude estimate. 

6.6. Dynamical interactions that postulate a high density of 
stellar black holes 

A dense cluster of stellar-mass BHs has been invoked as 
a potential solution to a number of other problems of col- 
lisional dynamics at the Galactic center Examples include: 
(1) removing stellar envelopes via physical collisions between 
BHs and stars (e.g. Dale et al . 2009); (2) randomizing the 
orbit s of young stars v ia gravitational scattering off of BHs 
(e.g. lPeretsetal.ll2009h ; (3) production of HVSs through en- 
counters with BHs (e.g. l O'Leary &Loebll2008l) . Typically 
«bh ^ '"^ is assumed, as in the relaxed, mass-segregated mod- 
els, implying a mass in BHs of ^ IO^Mq within one 0.1 pc. 
The 10—100 times lower BH densities found in some of the 
evolutionary models presented here would imply correspond- 
ingly lower rates of interaction. 

We note that alternative mechanisms exist for solving many 
of these outstanding problems. For instance, a single IMBH 
can randomize t he orbits of young s tars even more efficiently 
than a BH cusp jMerritt et al.ll2009l) . 

6.7. Cores and nuclear star clusters 

As shown above (§4), the relaxation time at the Galactic 
center is short enough that a parsec-scale core will shrink 
appreciably over the course of 10 Gyr. What is the connec- 
tion between such a hypothesized initial core, and the cores 
that are observe d in spheroids brighter than Mb ~ —19.5 
(iCdte et al.ll20'07l) ? 

The classification of spheroids into cored or coreless fam- 
ilies is based on data with an angular resolution of ^ 0.1", 
corresponding to a linear size of ^ 1 pc at a distance of the 
Virgo galaxy cluster. Neither the current core at the cen- 
ter of the Milky Way, nor the larger initial core postulated 
here, would be easily discernable at this distance. Instead, 
the Milky Way would likely be classified as a galaxy with a 
nuclear star cluster (NSC): its luminosity profile is relatively 
flat outside ~ 10 pc and rises steeply inside (.Launhardt et al] 
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l200llSchodeTetani2008tlGraham & Spitleril2009t) . The ob- 
served core sits atop that star cluster. NGC 205 also has both a 
nuclear star cluster and a core; the core radius is ~ 0.12" ss 0.5 
pc (|Merritt 2009), making it very similar to the Milky Way 
core (although it is not clear that NGC 205 contains a SMBH; 
Valluri et al. 2005). 

Many other galaxies with NSCs could also contain un- 
detected cores. Bright young stars, like those at the cen- 
ter of the Milky Way, M 31, and other galaxies with NSCs 
(Ivan der Marel et ani2007h . would tend to mask the existence 
of a core in the old population, as indeed they did until very 
recently at the Galactic center. 

While most NSCs are too small for their internal structure to 
be resolved, the half-mass relaxation time fnuc can be reliably 
estimated for many, and its mean dependence on host galaxy 
(not NSC) luminosity is 

logio(fnuc/yr) -9.38 -0.434 (Mb +16) (81) 

jMerrittl [20091) where Mb is the absolute blue magnitude of 
the bulge component. If we assume that galaxies with NSCs 
also contain SMBHs, and that the relaxation time at rinfi is 
no greater than its value at the NSC half-light radius, then 
relaxation times drop below 10 Gyr at Mb sa — 17, slightly 
fainter than t he estimated luminosity of the Milky Way bulge, 
Mb = -17.6 (i Marconi & Hun3l2003i) . 

These arguments suggest that cores comparable in size 
to SMBH influence radii might exist in other galaxies with 
NSCs. Relaxation times in the NSCs are expected to be short 
enough that such cores could shrink appreciably in 10 Gyr. 

7. CONCLUSIONS 

1. The distribution of old stars at the Galactic center ex- 
hibits a low-density core of radius ^ 0.5 pc. The deprojected 
central density is poorly constrained but is consistent with 
zero. 

2. Assuming that the old stars trace the mass in the inner 
parsec, the two-body relaxation time (for Solar-mass stars) 
is nowhere shorter than ^ 5 Gyr. The relaxation time at the 
influence radius of SgrA*, ring w 2.5 pc, is robustly estimated 



to be 20 - 30 Gyr. 

3. Reproducing the observed distribution of old stars with a 
steady-state distribution function requires a strongly truncated 
phase-space density at low energies and/or low angular mo- 
menta. If the stellar density increases more slowly than r^*^^ 
toward SgrA*, the velocity distribution must be anisotropic in 
the inner parsec, with a deficit of eccentric orbits ("anisotropic 
core"). Otherwise the distribution function can be isotropic 
("isotropic core"). 

4. Anisotropic core models evolve toward isotropy on a 
Gyr time scale. In the process, the core radius decreases only 
slightly. The observed (small) degree of anisotropy at the 
Galactic center is consistent with such models at both early 
and late times. 

5. On a longer time scale, gravitational encounters produce 
changes in stellar orbital energies, causing a pre-existing core 
to shrink. Initial core radii in the range 1 — 1 .5 pc evolve, after 
10 Gyr, to cores of the currently observed size. 

6. The dynamical friction force acting on an inspiralling 
massive body falls essentially to zero at roughly 1/2 the stellar 
core radius. This results in an accumulation of lOM© black 
holes in a shell just inside the stellar core. Orbital decay of 
an intermediate mass black hole would also be expected to 
stall at this radius, rather than the much smaller stalhng radius 
expected in a dense stellar cusp. 

7. The expected density of lOM© black holes in the inner 
parsec depends sensitively on their initial distribution and on 
the elapsed time, but may be substantially lower than in mod- 
els that assume the absence of a stellar core. Rates of gravi- 
tational wave driven inspirals of stellar-mass black holes are 
1-2 orders of magnitude lower than predicted by steady-state 
models with a mass-segregated density cusp. 
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APPENDIX 

A. Orbital distributions in a core around a SMBH 



In the context of "isotropic core" models, we derive the distribution of orbital eccentricities that one would measure at a point 
inside the core and near to the SMBH. 

We assume that r <C r/, <C nnfl, where r is the point of observation, rj, defines the truncation energy through Et, — i^ifb), and 
rinfl is the SMBH influence radius. 

Assuming a power law in space density outside the core, the phase space density is given by equation ( l50b , 

m=fo\Ery\ E>Et 

= 0, E<Eb (1) 



where we have assumed that the gravitational potential is 



i.e. that r < rinfl. 
The velocity space volume element is 



(2) 



d^v = 4nv^dvsmQdQ (3) 
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where 9 is the angle between v and r and < 9 < 7t/2. Using 



^-r + ^M-^' ,4a, 

J = rv sin 9 = ^1 -e^-^^ , (4b) 

\/-2E 



this becomes 



~ dade (5) 

2r^Vr a 

where vv = v'cos9. The distribution of eccentricities at r is then given by the integration over a, or 

dN f'i(l-e) da , , 

— ocW — /fl. (6) 

de avr 

The lower limit is the semi-major axis of an orbit of energy Ey. The upper limit on a corresponds to an orbit of eccentricity e 
with periapse at r. The integral is zero unless 

r->Y (7) 

\—e 2 

or 

e>l-2-. (8) 
rb 



Writing 

and /(fl) a^l^^"^, this becomes 



V? 



GM. 



2---(l-.^)- 
a r 



(9) 



^/Af rrl(i-e) a^l2-1 J. 

=NQe da , e>l-2— . (10) 

de Jrt/l /2-l-(l-e2)£ ri, 



Defining the new variable x ~ {I —e)j this takes on the simpler form: 



2r 

l>e>l . (11) 

rh 

Figure [19] shows N{e), normalized to unit total number, at radii of 

0.3,0.1,0.03,0.01 (12) 

times rf, assuming y = 1.8. Also shown are the mean eccentricities at each radius. N{e) approximates a delta function at e = 1 
for r ^ rh, since the only orbits that reach into these small radii must be very eccentric. 
Proceeding as before, the distribution of orbital semi-major axes at r is 



1/2-Y 

, a > rb/2. (13) 



where e ~ I and r rt, have also been assumed. The mean value of a is 

{a) _ l y-3/2 
rb 25/2- J 

andforY== 1.8, 



(14) 



(a) 

— =0.214. (15) 
rb 

Setting rb ~ 0.5 pc, the expected orbit should have a semi-major axis of ^ 0.1 pc and an eccentricity near one. 
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